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I. 

INTRODUCTION 

This  report  summarizes  research  carried  out  under 
the  Post-Doctoral  Program  in  Seismology  during  the  period 
1  July  1971  to  30  June  1972  and  earlier  work  that  was  sub¬ 
mitted  for  publication  during  this  period.  Staff  members 
of  the  Department  of  Earth  and  Planetary  Sciences  and 
Lincoln  Laboratory  participated  in  this  program. 

The  work  described  here  is  divided  into  three 
categories:  (1)  properties  of  seismic  sources,  (2)  studies 

of  the  dynamical  behavior  of  the  crust  and  upper  mantle,  and 
(3)  observations  of  crustal  structure. 

The  focal  depths  of  the  three  largest  earthquakes  in 
the  Parkfield,  California,  sequence  were  determined  using  first 
and  second  P  arrivals  at  Berkeley,  California.  Data  were 
insufficient  to  determine  whether  the  second  arrival  corre¬ 
sponds  to  refraction  along  a  discontinuity  within  the  crust 
or  to  a  complicated  source-time  function.  With  either 
interpretation,  however,  the  focal  depths  for  the  three 
earthquakes  apparently  increased  with  time  through  the  sequence. 

Numerical  calculations  have  been  made  of  spherical  wave 
propagation  due  to  explosions  in  a  number  of  rock  types.  The 
effect  on  the  seismic  source-time  function  of  small  nonlinear 
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irreversibility  is  examined  for  propagation  at  distances 
where  the  response  is  commonly  regarded  as  elastic.  Amplitude 
of  long-period  components  is  found  to  increase  due  to  non¬ 
linear  coupling  between  frequencies. 

A  numerical  simulation  of  upper  mantle  convection  has 
been  constructed  using  a  two-dimensional  time -dependent 
model  that  allows  large  viscosity  variations.  The  method 
has  been  applied  to  sea- floor  spreading  with  the  objective 
of  examining  the  driving  mechanism  of  mid-ocean  ridges. 
Counterflow  below  the  plates  is  confined  to  depths  less  than 
340  km.  It  is  found  that  for  a  spreading  velocity  of  1.2 
cm/yr,  the  ridge  can  produce  compressive  stress  in  the 
lithosphere  out  to  a  distance  of  1600  km.  For  a  spreading 
rate  of  6  cm/yr,  however,  this  model  is  clearly  excluded 
for  it  requires  an  excessively  large  stress  in  the  lithosphere. 
Therefore,  upwelling  material  must  cross  the  seismic  dis¬ 
continuity  at  the  400  km  depth. 

A  numerical  method  for  more  general  rheological  models 
has  also  been  developed. 

The  crustal  structure  beneath  LASA  has  been  investigated 
using  the  spectral  ratio  of  vertical-to-horizontal  displace¬ 
ments  of  long-period  P  waves.  The  period  of  the  lowest 
frequency  peak  in  the  spectral  ratio,  related  to  the  vertical 
P  wave  travel-time  in  the  crust,  is  quite  variable  over  the 


3 


area  spanned  by  the  array.  Variations  in  crustal  thickness 
of  about  7  km  are  implied.  The  Moho  generally  shoals  from 
the  northeast  to  the  southwest  across  the  array  and  exhibits 
a  synclinal  structure  with  the  axis  plunging  toward  the  north¬ 
east  in  the  southwest  quadrant  of  LASA. 

Some  of  the  work  discussed  herein  has  already  been  pub¬ 
lished.  The  summary  for  such  work  is  given  below  in  abstract 
form.  The  complete  results  may  be  found  by  consulting  the 
appropriate  reference  given  in  Section  V.  of  this  report. 
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II.  SEISMIC  SOURCES 

11. 1  Focal  Depths  of  the  1966  Parkfield,  California, 
Earthquakes  by  William  H.  Bakun  (Abstract) 

Differences  in  arrival  times  of  seismic  phases  at 
Berkeley,  California  (BRK) ,  A  -  270  km,  for  the  0408 
UT  June  28  (M  =  5.1),  0426  UT  June  28  (M  =  5.5),  and 
1953  UT  June  29  (M  =  5.0)  1966  Parkfield,  California, 
earthquakes  imply  that  focal  depths  for  these  three  largest 
events  of  the  1966  Parkfield  sequence  increased  with  time 
through  the  sequence.  The  data  available  are  not 
sufficient  to  determine  whether  the  observed  secondary 
arrivals  at  BRK  result  from  a  slower  propagation  path  or 
are  part  of  a  complicated  source-time  function. 

11. 2  Propagation  of  Underground  Explosion  Waves  in  the 
Nearly-Elastic  Range  by  D.H.  Andrews  and  Seymour 
Shlien 

Abstract 

The  effect  of  small  anelasticity  on  the  propagation  of 
spherically  symmetrical  waves  is  examined  with  numerical 
calculations.  The  type  of  anelasticity  used  is  static  hyster¬ 
esis  with  nominal  Q  value  of  100.  Linear  theories  do  not  apply 
in  this  case.  Fourier  components  of  different  frequencies 
do  not  propagate  independently,  but  energy  is  transferred  to 
lower  frequencies.  The  extent  to  which  reduced  displacement 
potential  is  not  invariant  with  respect  to  radius  is  examined 
for  explosions  in  granite,  salt,  and  shale. 
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In  this  paper  nonlinear  effects  on  propagation  of  pulses 
from  underground  explosions  are  examined  in  the  nearly  elastic 
region.  The  calculational  method  used  is  similar  to  numerical 
methods  used  in  AEC  laboratories  (Rodean,  1971)  for  propagation 
in  the  close-in  region,  where  pulse  widths  change  from 
microseconds  to  milliseconds.  We  are  concerned  with  propaga¬ 
tion  beyond  the  radius  at  which  the  material  response  is 
commonly  considered  to  be  elastic,  and  will  consider  the  effect 
of  small,  but  permanent,  irreversible  strains. 

The  mechanism  of  attenuation  in  crustal  rocks  has  been 
reviewed  by  Knopoff  (1964)  and  by  Gordon  and  Davis  (1968) . 

They  find  that  in  polycrystalline  samples  the  stress-strain 
path  is  not  reversible,  but  has  a  small  hysteresis  that  is 
nearly  independent  of  strain  rate.  This  effect  is  not 
plasticity  in  the  usual  sense,  for  the  loading  and  unloading 
slopes  are  only  slightly  different.  The  effect  probably 
arises  from  friction  at  grain  boundaries. 

Knopoff  and  MacDonald  (1958)  and  Carpenter  (1966)  have 
developed  linear  viscoelastic  models  for  which  Q  is  constant 
over  a  broad  frequency  band.  These  models  are  rather  com¬ 
plicated.  Their  work  was  motivated  by  trying  to  maintain  the 
principle  of  superposition,  so  that  wave  propagation  could 
be  treated  analytically.  However,  it  is  possible  that  a  con¬ 
stant  Q  could  be  due  to  static  hysteresis.  This  is  the 
most  plausible  type  of  hysteresis  if  the  energy  loss  is  due 
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to  friction  at  grain  boundaries. 

Johrison  (1955)  measured  the  response  of  an  assembly 
of  elastically  loaded  bodies  with  some  sliding  at  their 
interfaces,  and  found  a  stress-strain  hysteresis  loop  that 
was  the  same  even  when  frequency  approached  zero. 

In  the  case  of  static  hysteresis,  the  principle  of 
superposition  does  not  hold,  and  different  frequency  com¬ 
ponents  will  not  propagate  independently.  For  waves  from 
underground  explosions,  energy  transfer  from  higher  to  lower 
frequencies  could  be  important. 

Since  we  wish  to  calculate  wave  propagation  with  non¬ 
linear  material  properties,  a  numerical  method  must  be  used. 

We  adopt  the  finite  difference  equations  of  Wilkins  (1964) 
in  which  the  first  principle  equations  for  continuum  motion 
are  approximated  directly.  These  equations  for  the  case  of 
spherical  symmetry  are  given  in  the  Appendix. 

In  the  material  model  used,  the  shear  response  is  ir¬ 
reversible,  as  shown  in  Figure  1.  The  hysteresis  is  independ¬ 
ent  of  strain  rate.  When  shear  stress  and  shear  strain  rate 
have  the  same  sign  (loading) ,  a  constant  shear  modulus  Uq  is 
used.  When  the  stress  and  the  strain  rate  have  opposite 
signs  (unloading) ,  stress  increments  are  related  to  strain 
increments  using  a  variable  shear  modulus, 

u  .  „  +  (1) 
max 

where  s_  is  the  maximum  shear  stress  reached  at  the  par tic- 
max 
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ular  locality  in  the  previous  loading*  Therefore,  the 
initial  unloading  slope  is  larger  than  the  loading  slope, 
and  then  decreases  smoothly  to  the  same  value.  Behavior  is 
similar  for  loading  and  unloading  in  the  opposite  direction. 

In  this  model  the  residual  strain  when  stress  is  reduced 
to  zero  is  proportional  to  the  maximum  stress  reached.  The 
relative  energy  loss  is  independent  of  amplitude  and  of 
frequency.  The  principle  of  superposition  does  not  apply. 

We  will  calculate  spherical  wave  propagation  for  the 
explosions,  Gasbuggy,  29  kilotons  detonated  in  shale,  Hardhat, 
5  kilotons  in  granite,  and  Salmon,  5  kilotons  in  salt.  In 
'■^iilT^t'iTj^the  Vaveform  is  specified  at  a  radius  beyond  the 
strongly  inelastic  region,  and  propagation  to  greater 
distances  is  calculated.  It  is  not  reasonable  to  assume  that 
the  material  is  the  same  at  these  distances  as  at  the  shot 
point.  The  Pictured  Cliffs  Sandstone  formation  in  which  a 
Gasbuggy  gauge  was  located  is  chosen  as  a  typical  near  surface 
crustal  rock.  Its  properties  are  used  in  all  calculations. 
Density  is  2.51  gm/cm  ,  P  wave  velocity  is  4.3  km/sec,  and 
Poisson's  ratio  is  0.3  (Perret,  1969).  In  the  inelastic 
calculations,  the  shear  modulus  was  calculated  from  equation 
(1)  with  u0  *  0.1327  megabar  and  Ay  -  0.0328  megabar.  This 
corresponds  to  a  relative  energy  loss  in  shear  of  16% ,  or  in 
uniaxial  strain  of  6%,  giving  a  nominal  Q  for  P  waves  of 
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The  Gasbuggy  waveform  at  468  meters,  with  a  peak  velocity 
of  2,4  meter/sec,  peak  displacement  of  6.8  cm,  and  final 
displacement  of  2.6  cm,  (Perret,  1969)  can  be  matcned  roughly 
(in  the  elastic  case)  by  a  stress  pulse  of  the  form 

Ex  *  P0  +  Pie”t/T  (2) 

For  the  Gasbuggy  calculations  units  of  length  and  time 
are  scaled  down  by  the  factor  (5/29)^^,  so  that  all  calcul¬ 
ations  are  for  5  kilotons.  values  of  the  parameters  in 
equation  (2)  for  Gasbuggy  scaled  to  5  kilotons  are  given  in 
Table  I.  The  Salmon  waveform  at  625  meters,  with  a  peak  veloc¬ 
ity  of  2  meter/sec,  peak  displacement  of  2.8  cm,  and  final 
displacement  of  0.85  cm  (Patterson,  1964)  is  also  fit  by  a 
stress  pulse  of  the  form  given  in  equation  (2)  with  parameters 
listed  in  Table  I. 

The  Hardhat  calculation  is  not  based  on  a  measured 
waveform.  There  was  a  significant  amount  of  tectonic  strain 
release  in  the  Hardhat  event  (Toksoz  and  Kehrer,  1971) ,  which 
could  have  contributed  to  the  measured  displacement.  Our 
calculation  is  based  on  the  stress  pulse  calculated  by  Cherry 
(Rodean,  1971) ,  which  can  be  fit  closely  by  equation  (2) 
with  parameters  listed  in  Table  1. 

In  order  to  facilitate  comparison  of  attenuation  over 
equal  distances,  it  was  decided  to  do  all  the  numerical 
calculations  with  the  input  waveform  imposed  at  the  same 
radius,  400  meters.  The  stress  profiles  specified  in  Table  1 
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are  transformed  by  varying  pQ  inversely  with  the  cube  of 
radius  to  get  the  same  final  displacement  and  varying  p^ 
inversely  with  radius  to  get  the  same  peak  velocity  as  a 
function  of  radius  in  the  elastic  case.  The  transformed 
parameters  are  listed  in  Table  II. 

Although  this  transformation  of  the  stress  pulse  keeps 
peak  velocity  and  final  displacement  invariant,  it  changes 
the  spectrum  of  the  wave  generated.  In  the  elastic  case  the 
period  of  the  predominant  low  frequency  component  is  proportion¬ 
al  to  the  radius  at  which  the  stress  pulse  of  equation  (2) 
is  applied  (Sharpe,  1942) .  We  have  not  taken  care  to  establish 
realistically  the  radii  at  which  the  three  different  ex¬ 
plosion  pulses  become  nearly  elastic,  but  we  will  investigate 
the  effect  of  nonlinear  material  response  on  three  waveforms 
with  varying  ratio  of  impulse  to  step  function. 

For  each  explosion  an  elastic  and  an  inelastic  calculation 
was  done.  In  the  elastic  calculations,  the  stress  pulse 
specified  in  Table  II  was  applied  at  400  meters.  The  calculated 
method  was  checked  by  verifying  that  the  reduced  displace¬ 
ment  potentials  calculated  at  different  radii  were  the  same. 
Velocity  as  a  function  of  time  at  meters  was  saved  from 
each  of  the  elastic  calculations  and  was  used  as  a  boundary 
condition  for  the  inelastic  calculations.  The  finite  dif¬ 
ference  zone  size  in  all  calculations  was  5  meters. 
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Velocity  pulses  calculated  for  Gasbuggy  at  a  radius  of 
1.4  km  are  compared  In  Figure  2.  The  waveform  for  the 
inelastic  case  looks  similar  to  the  elastic  waveform  but  has 
smaller  amplitude.  An  important  feature  not  evident  in 
this  figure  is  that  the  final  displacement  is  larger  in  the 
inelastic  case. 

Reduced  displacement  potential  is  calculated  from  the 
equation 

Y(t)  =  cr  Jfc  d(t')ec(t“t,,/rdt'  (3) 

0 

where  d  is  displacement.  Potentials  for  the  three  explosions 
in  the  case  of  elastic  propagation  beyond  400  meters  are 
shown  as  solid  curves  in  Figure  3. 

In  the  inelastic  case  the  reduced  displacement  potential 
has  no  meaning,  although  the  integral  in  equation  (3)  can 
be  performed.  Indeed,  this  is  what  is  done  when  a  "reduced 
displacement  potential"  is  derived  from  measurements ,  for 
there  is  no  assurance  that  the  material  response  is  strictly 
elastic.  In  the  inelastic  calculations  y  calculated  from 
equation  (3)  at  400  meters  is  the  same  as  the  elastic  case, 
since  displacement  is  prescribed  to  be  the  same  there.  At 
greater  distances,  however,  ¥  is  not  the  same.  Dashed  curves 
in  Figure  3  show  y  for  the  three  explosions  at  1.4  km.  The 
difference  between  the  solid  and  dashed  curves  is  the  dif¬ 
ference  that  might  be  expected  in  "potentials”  derived  from 
measurements  1  km  apart  for  spherical  wave  propagation  in 
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rock  with  Q  =»  100.  The  initial  peak  of  ¥  is  smaller 
due  to  attenuation  of  the  propagating  pulse,  but  then  ¥ 
rises  to  larger  values  because  of  the  larger  final  disDlace- 
ment. 

These  inelastic  "reduced  displacement  potentials" 
can  not  be  rigorously  related  to  more  distant  seismic  waves, 
but  some  qualitative  conclusions  can  be  drawn.  If  the 
material  had  been  changed  from  inelastic  to  elastic  at  1.4  km, 
the  reduced  displacement  potential  calculated  there  would 
have  been  rigorously  meaningful.  In  this  case  the  initial 
peak  of  ¥  would  not  change  much,  since  it  is  determined  by 
attenuation  within  1.4  km.  Later  values  of  ¥  would  be  smaller 
since  the  final  displacement  is  determined  by  yielding  both 
within  and  beyond  1.4  km.  However,  these  later  value  of  ¥ 
would  still  be  larger  than  in  the  all  elastic  case. 

The  functions  shown  in  Figure  3  are  analyzed  here  as  if 
they  were  truly  seismic  source-time  functions.  Particle 
velocity  of  surface  waves  and  head  waves  is  proportional  to 
the  first  derivative  of  ¥.  Fourier  amplitude  spectra  of 
d¥/dt  are  shown  in  Figure  4  for  the  elastic  and  inelastic 
cases  at  1.4  km  for  the  three  explosions.  The  location  of 
the  spectral  peak  near  2  hertz  is  determined  by  the  choice 
of  400  meters  as  the  cavity  radius.  At  this  frequency  and 
above  the  amplitude  in  the  inelastic  cases  is  reduced  due 
to  attenuation  through  1  km  of  irreversible  material.  The 
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The  DC  component  is  larger  due  to  the  larger  final  displacement. 

i 

This  would  not  be  the  case  with  a  linear  viscoelastic  atten¬ 
uation  model.  Since  the  principle  of  superposition  does  not 
apply  in  this  model ,  different  components  do  not  propagate 
independently.  Energy  can  be  transferred  from  higher  to  lowec 
frequencies.  In  the  spectra  shown  it  is  evident  that  the 
increased  energy  in  the  DC  component  is  taken  from  the  1  hertz 
component.  This  minimum  might  change  as  the  wave  propagates 
farther. 

The  particle  velocity  of  distant  body  waves  is  propor¬ 
tional  to  the  second  derivative  of  t.  Therefore,  body  wave 
spectra  are  obtained  by  multiplying  surface  wave  spectra  by 
frequency.  Therefore,  the  inelastic  effects  at  low  frequency 
will  be  less  dramatic  for  body  waves  than  for  surface  waves. 

The  quality  factor,  Q,  of  100  used  here  is  reasonable 
for  near-surface  crustal  rocks.  Also,  static  hysteria  is 
consistent  with  laboratory  measurements.  The  fact  that 
static  hysteresis  has  not  been  considered  before  in  seismological 
theory  is  related  more  to  difficulty  of  analysis  than  to 
indications  of  the  data.  Therefore,  the  inelastic  model  used 
here  is  a  reasonable  possibility.  The  seismic  source-time 
function  for  an  underground  explosion  can  be  significantly 
different  from  that  derived  at  the  commonly  accepted  elastic 
radius. 
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TABLE  1 


Radius 

P0 

pl 

T 

Gasbuggy  (scaled) 

260  m 

29.7  bar 

259  bar 

15.7  msec 

Hardhat 

400 

11 

144 

10 

Salmon 

625 

7.2 

215 

14.1 
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TABLE  II 


Radius 

P0 

Pi 

T 

Gasbuggy  (scaled 

400  m 

8.15  bar 

168  bar 

15.7  msec 

Hardhat 

400 

11 

144 

10 

Salmon 

400 

27.5 

336 

14.1 

Figure  Captions 


Figure  1: 


Figure  2: 


Figure  3: 


Figure  4: 


A  typical  path  for  shear  stress  versus 
shear  strain  for  the  inelastic  model. 

Particle  velocity  at  1.4  km  for  the 
Gasbuggy  calculations.  Solid  curve  is 
the  elastic  case  and  dashed  curve  is 
the  inelastic  case. 

Solid  curves  are  reduced  displacement 
potentials  for  the  three  explosions  in 
the  elastic  case.  Dashed  curves  are 
functions  derived  by  equation  (3)  from 
displacement  at  1.4  km  in  the  inelastic  case. 

Fourier  amplitude  spectra  of  surface  wave  source  time 
functions.  Solid  curves  are  the  elastic  case. 

Dashed  curves  are  found  by  taking  dashed 
curves  of  Fig.  3  as  reduced  displacement 
potentials. 
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APPENDIX 

A  Lagrangian  finite  difference  scheme  is  used.  Space 
intervals  are  indexed  by  the  subscript  j  and  time  steps  are 
indexed  by  the  superscript  n.  Grid  points  move  with  the 
material,  so  that  each  interval  between  points  always  con¬ 
tains  the  same  mass. 

The  equation  of  motion 

3u  .  1  3Ir  .  ,  rr  -  h 
Tt  p  Tr  p  r 

is  approximated  by  the  finite  difference  equation 

-n*r>;tV2  -'VS-1/21^  +  2B” 


where 


[»j  +  1/2  (l£l"  rj>  +  pj-l/2(5n*  rj-l' 


+  1/2  ~  +  1/2  +  (£r)i-l/2~(£8)  j-l/2_ 

1/2  p4  +  1/2  <rj  +  1  +  rj’  1/2  Pj-l/2(r”  +  rj-l’ 


Then  acceleration  is  integrated  to  update  velocity 
un+l/2  .  un-l/2  +  At 

and  the  velocity  is  integrated  to  update  the  radius  of  each 
Lagranian  point 


r!j+1  -  r"  +  At  u"+1/2 
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Updated  velocities  and  co-ordinates  are  used  to  evaluate 
strain  rates. 


n+1/2 

j+1/2 


$+1/2  -  <on+1/2  ♦  Jlijsw 


j+1/2 


r  j+1/2 


e' j+1/2 


Then  a  new  value  of  pressure  is  found 

Pn+1  .  Pn  .  k  (V}n+1/2  At 
P j+1/2  P j+1/2  k  v  j+1/2  At 

where  k  is  bulk  modulus.  Updated  stress  deviators  are 


given  by 


(sr} j+1/2  =  (s 


>n  +  2u  Ff f  \n+l/2  1  v"+1/21 

r  j+1/2  +  2y|(er)  j+1/2  '  7  (V>+1/2  J 


/_  >n+l  1.  .n+1 

j+1/2  ■  7'sr' j+1/2 


Irreversible  behavior  of  the  stress  deviator  s  is 

r 

obtained  if  the  shear  modulus  y  used  in  the  incremental 
equation  above  is  not  constant#  but  depends  on  sr  and  on 
the  direction  of  change  of  sr- 
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Total  stress  components  are 

<Er> ?i/a  *  <-p  +  *r>jtv2 

(e8>  "ti/2 =  (-p  +  vSU 


After  this  sequence  of  equations  is  executed  for  each 
value  of  j,  the  time  step  index  n  can  be  incremented  by  one 
and  the  cycle  can  be  repeated.  All  equations  are  explicit. 
The  stability  requirement:  is 

it  <  ^ 

c 

where  c  is  sound  speed. 
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III.  DYNAMICS  OF  THE  CRUST  AND  UPPER  MANTLE 

111.1  Numerical  Simulation  of  Sea-Floor  Spreading 
by  D.J.  Andrews  (Abstract) 

Upper  mantle  convection,  including  lithospheric  plate 
motion,  is  simulated  numerically  using  a  two-dimensional 
time-dependent  method  that  allows  large  viscosity  variations. 
The  numerical  operator  developed  for  viscous  flow  is  in  self- 
adjoint  form,  so  that  the  conjugate  gradient  iteration  method 
may  be  used.  Convergence  is  faster  than  with  relaxation 
methods.  The  method  is  applied  to  sea-floor  spreading  with 
the  objective  of  examining  the  driving  mechanism  of  mid-ocean 
ridges.  In  accordance  with  this  objective,  deep  convection 
is  suppressed  in  the  model.  Counter flow  below  the  plates 
is  confined  to  depths  less  than  340  km.  The  model  is  fit  to 
observed  topography  at  four  different  ridge  locations.  It 
is  found  that  for  a  spreading  velocity  of  1.2  cm/yr,  the 
ridge  can  produce  compressive  stress  in  the  lithosphere  out 
to  a  distance  of  1600  km.  For  a  spreading  velocity  of  6  cm/yr 
this  model  is  clearly  excluded,  for  it  requires  excessively 
large  stress  in  the  lithosphere.  Therefore  upwelling  material 
must  cross  the  seismic  discontinuity  at  400  km  depth. 

111. 2  A  Numerical  Method  for  Creep  Deformation  of  Solids 
by  D.J.  Andrews 
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This  note  is  concerned  with  extension  of  the  procedure 
described  by  Andrews  and  Hancock  [1]  to  time-dependent  problems 
in  which  motion  is;  slow  enough  that  inertial,  forces  are 
negligible.  Time  dependence  of  the  solution  may  arise  from 
stress  relaxation  of.  the  material  and  from  time-dependent 
boundary  conditions.  The  procedure  involve.*  advancing 
through  time  in  finite  steps  and  iterating  io  acheive  stress 
equilibrium  at  each  step  in  time.  It  is  discussed  in  terms 
of  a  Maxwellian  vi scoelastic  model  that  can  be  generalized 
to  nonlinear  case?;. 

Stress  components  0^  are  decomposed  into  stress  de¬ 
viators  and  pressure 


°ij  =  sij  -  p  «ij 

Pressure  is  uniquely  determined  by  volume,  but  stress  de¬ 
viator  components  obey  a  stress  relaxation  law,  which  in 
differential  form  is 


!2ji  'dc. 


-  s .  ■ 
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dt/T 


where  e. .  is  strain  deviator,  y  is  shear  modulus,  and 
t  is  relaxation  lime.  If  x  is  const; nt  this  is  linear 
Maxwellian  viscoelasticity,  in  non line;  r  cases  x  is  a 
function  of  stress..  To  have  a  properly  covariant  descrip¬ 
tion,  x  should  be  expressed  as  a  function  of  stress  in- 
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variants. 

The  finite  difference  equation  used  to  advance  stress 
in  one  zone  from  time  step  n  to  step  n  +  1  is  derived 
as  follows.  Let  (*i;j)n+1^2  be  the  strain  deviator  rate, 
found  from  velocities,  for  that  zone  for  advancing  from 
time  step  n  to  a  particular  iteration  at  time  step  n  +  1. 

A  finite  difference  analogue  of  the  differential  -equation 
above  is 

<s..)n+1  =  (s.  .)"  +  2u(4.j)n+1',2  it 

-  l/2[(s1;j)n+1  +  <s..)n]  |t 

This  may  be  rearranged  to  get  an  explicit  equation 

<sij)n+1  =  +  2p (&i j ) n+1/2  4t 

-  l/2(s..)n  J/U  +  1/2 

This  equation  is  stable  for  all  values  of  At  and  is 
accurate  to  second  order  in  At/t 

An  iteration  must  be  performed  to  converge  to  stress 
equilibrium  at  time  step  n  +  1.  Solution  of  the  above 

equation  for  all  zones  constitutes  step  4  of  the  iteration 
outlined  in  [1] . 
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To  proceed  through  the  next  iteration  at  time  step 
n  +  1  ,  stress  components  just  calculated  are  used  to  find 
the  unbalanced  force  on  each  grid  point  (step  1  of  the  itera 
tion) .  Then,  grid  points  are  displaced  in  the  direction  of 
this  force  to  go  from  positions  in  the  previous  iteration 
at  time  step  n  +  1  to  positions  in  the  current  iteration 
at  time  step  n  +  1  (step  2) .  The  velocities  of  grid 
points  from  time  step  n  to  the  current  iteration  in  time 
step  n  +  1  ,  are  found.  Prom  these  velocities  strain 
rates  are  found  (step  3) ,  and  then  the  stress  calculation 
may  be  repeated. 

In  the  case  of  nonlinear  stress  relaxation,  the  re¬ 
laxation  time  t  should  be-  evaluated  from  invariants  of 
the  stress  tensor  averaged  at  the  new  and  old  times.  In 
this  average  one  may  use  stress  at  time  step  n  and  stress 
from  the  previous  iteration  at  time  step  n  +  1. 

This  procedure  has  been  used  in  a  problem  with  a  cubic 
creep  law.  The  iteration  behaved  in  a  reasonable  way. 

To  check  the  accuracy  of  the  method  a  problem  was  done 

4 

with  a  linear  viscoelastic  material  in  an  infinite  half 
space,  with  a  pressure  applied  to  the  surface.  The  x-axis 
extends  into  the  medium  and  the  surface  is  at  x  =  0.  The 
material  has  been  at  rest  with  no  pressure  on  the  surface 
at  all  times  ;up  to  t  =  0.  At  t  =  o  the  pressure 

p  =  P  cos  ay 
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is  suddenly  applied  to  the  surface  and  is  held  constant 
thereafter.  The  analytic  solution  is  found  by  applying 
Bland's  correspondence  principle  [2]  to  the  elastic 
solution  [3] .  The  components  of  displacement  are 

u  =  e_aX  cos  a  y  [-  £  U  -  a)  e“at/T 

+  £  +  1  +  ax 

+  (1  +  ax)  t/t] 

and 

v  “  ■  2uo  e"OIt  sin  *rl-f  ll-al  e-at/l 

+  k  '  ox 
-  ax  t/x] 

where  k  is  bulk  modulus  and 
a”1  “  1  +  y/(3k) 

At  t  =0  these  expressions  are  the  solution  for  the 
elastic  case  [3].  Note  that  in  the  elastic  case  the  hori¬ 
zontal  displacement  reverses  direction  at  a  depth 

ax  =  1  -  2v  , 

while  in  the  viscoelastic  solution  the  horizontal  velocity 
at  late  times  is  in  the  same  direction  at  all  depths. 
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In  the  numerical  test  case  we  will  choose  t  =  1, 
<**1,  2y  =  l,  v  =  0.2.  The  numerical  method  is  valid 
for  large  displacements ,  but  the  analytic  solution  holds 
only  for  small  displacements .  To  keep  displacement  small 
Ve  '  *ioose  P  -  10  *.  Displacements  are  multiplied  by  10* 
in  the  figures. 

In  the  finite  difference  calculation  8  zones  are 
used  in  a  half  wavelength  of  the  pressure  variation,  and 
the  region  considered  is  12  zones  deep.  Zones  are  approx¬ 
imately  square.  The  pressure  was  suddenly  applied  at  time 
zero,  and  the  calculation  proceeded  for  one  relaxation  time. 
The  time  step  used  was  t/10.  In  each  time  step  200 

» 

iterations  were  performed  to  approach  stress  equilibrium. 

The  number  of  iterations  required  for  long  wavelength  com¬ 
ponents  to  converge  increases  when  finer  zoning  is  used. 

-It  is  proportional  to  the  square  of  the  number  of  zones  in 
one  dimension. 

Displacements  calculated  afLer  the  first  time  step  are 
shown  in  Figure  1.  This  is  approximately  the  elastic  so¬ 
lution  expected  for  instantaneous  displacement.  Variation 
of  each  component  of  displacement  in  the  direction  parallel 
to  the  surface  is  sinusoidal,  as  it  should  be,  within  1 
percent.  Time  dependence  of  the  x-component  of  displacement 
at  y  =  0  is  shown  for  four  different  depths  in  Figure  2. 
Symbols  show  calculated  values  and  the  solid  curves  are  the 
analytic  solution.  In  the  first  time  step  errors  are  6 


percent  of  the  maximum  displacement.  This  error  could  have 
been  reduced  by  using  a  larger  number  of  iterations.  The 
deviation  from  stress  equilibrium  is  partly  corrected  in 
the  next  time  step,  where  errors  are  less  than  2  percent 
of  maximum  displacement.  Time  dependence  of  the  y-component 
of  displacement  at  y  =  tt/2  is  shown  in  Figure  3. 
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Figure  1 


Figure  2 


Figure  3 . 


Displacement  field  in  the  demonstration  problem 
after  the  first  step  in  time. 


Vertical  displacement  at  four  different  depths 
as  a  function  of  time.  Symbols  are  calculated 
values.  Solid  curves  are  the  analytic  solution. 

Horizontal  displacement  at  four  different  depths 
or  a  function  of  time.  Symbols  are  calculated 
values,  solid  cuives  are  the  analytic  solution. 
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VERTICAL  DISPLACEMENT 
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IV.  STRUCTURE  OF  THE  CRUST 


IV. 1  Crustal  Structure  Beneath  LASA  from  Long-Period 
P-Wave  Spectra  by  William  H.  Bakun  (Abstract) 

Long-period  transfer-function  ratios  from  events  in 
South  America,  Fiji-Tonga,  and  Japan  recorded  at  the  LASA 
subarray  centers  are  interpreted  in  terms  of  Haskell- 
Thomson  theory.  The  transfer-function  ratio  data  provide 
>  three-dimensional  model  for  the  crustal  structure  beneath 
LASA.  The  proposed  structure  can  be  characterized  by  two 
trends:  (1)  crustal  thinning  from  the  northeast  to  the 

southwest  across  the  array  and  (2)  a  synclinal  structure 
in  the  southwest  quadrant  of  the  array  with  axis  plunging 
toward  the  northeast. 
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